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Abstract 

We determine the non-linear hydrodynamic response to geometrical fluctuations in heavy ion 
collisions using ideal and viscous hydrodynamics. This response is characterized with a set of non- 
linear response coefficients that determine, for example, the v§ that is produced by an t2 and an €3. 
We analyze how viscosity damps both the linear and non-linear response coefficients, and provide 
an analytical estimate that qualitatively explains most of the trends observed in more complete 
simulations. Subsequently, we use these non-linear response coefficients to determine the linear and 
non-linear contributions to v\ , and v 5 . For viscous hydrodynamics the non-linear contribution is 
dominant for V4, v§ and higher harmonics. For v\, the non- linear response constitutes an important 
~ 25% correction in mid-central collisions. The non-linear response is also analyzed as a function 
of transverse momentum for v±, i>4 and ^5. Finally, recent measurements of correlations between 
event-planes of different harmonic orders are discussed in the context of non-linear response. 
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I. INTRODUCTION 



The goal of the RHIC and the LHC heavy ion programs is to produce and to characterize 
the Quark Gluon Plasma (QGP), a prototype for non-abelian plasmas. One of the best ways 
to understand the transport properties of the experimentally produced plasma is through 
anisotropic flow [1-3]. In a heavy ion collision the nuclei pass through each other, and the 
resulting energy density in the transverse plane fluctuates in coordinate space from event to 
event. If the mean free path is short compared to the system size, the produced plasma will 
respond as a fluid to the pressure gradients and convert these coordinate space fluctuations 
to long range momentum space correlations between the produced particles. In the last two 
years it was gradually realized [4-6] that all of the long range momentum-space correlations 
known colloquially as the "ridge" and "the Mach cone" are manifestations of this collective 
flow [7, 8]. This realization gave rise to a large variety of flow observables which provide 
an unprecedented experimental check of the overall correctness of the hydrodynamic picture 
of heavy ion events [7, 9-11]. Further, different observables have different sensitivity to 
the shear viscosity of the plasma [12], and therefore a global analysis of flow can provide 
cross-correlated constraints on r]/s. 

One of the most direct measurements is the harmonic spectrum of the produced particles. 
The final state momentum spectrum for each event can be expanded in harmonics 



where P is the azimuthal angle of the produced particles and \l/ ra is the event plane angle 1 . 
The averaged square of these harmonics, i.e. ((v^)), can be measured experimentally by 
studying two particle correlations [1]. There is strong experimental and theoretical evidence 
that the harmonic coefficients, v% and v$, are to a good approximation linearly proportional 
to the deformations in the initial energy density in the transverse plane. For example, the 
experimental ratio / ((i>f)) closely follows the geometric deformations ((e§)) / ((e^)) as a 
function of centrality [7]. Event-by-event simulations with ideal hydrodynamics reproduce 
this trend, and show that the event plane angles ^2 and ^3 are strongly correlated with the 
angles of the initial deformations [13]. 

However, in an insightful paper Gardim et al [14] studied the correlation between higher 
harmonics, V4 and v$, and the initial spatial deformations within ideal hydrodynamics. This 
work explained and quantified the extent to which the higher harmonics such as V4 and 
t>5 arise predominantly from the non-linearities of the medium response. For example, for 
mid-central collisions the observed v$ is predominantly a result of the interactions between 
t>2 and V3. This work was motivated in part by previous event-by-event simulations by 
Heinz and Qiu [13] which showed that ^4 and ^5 are uncorrelated with the fourth and 
fifth harmonics of the spatial deformation. Based on the centrality dependence of this 
decorrelation, these authors anticipated (but did not quantify) the importance of V2-V3 mode- 
mixing in determining V5. 

The goal of this work is to systematically characterize the non-linear response of the 
medium. First, in Section II we introduce a set of non-linear response coefficients, and 

1 Following tradition, we have expanded the particle distribution in terms of cosines and phases *S/ n rather 
than cosines and sines. 




(1.1) 
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describe how these coefficients can be used in conjunction with a Glauber model to determine 
((f^)). The strongest non-linear response stems from the interactions between v-i and the 
other harmonics, and consequently a prominent response coefficient is W5(23)/ e 2e3 ) which 
determines the V5 produced by an elliptic and triangular deformation. In Section III B we 
determine these response coefficients using both ideal and viscous hydrodynamics, and study 
how the response depends on the shear viscosity. With these non-linear coefficients, together 
with the linear response, we make several predictions for Vi, v±, and v$ in ideal and viscous 
hydrodynamics in Section IV. Finally, in Section IV we also study the transverse momentum 
dependence of v%, t>4, and v$. 

In this work we will determine the harmonic spectrum by characterizing the quadratic 
response of the system to small deformations. Alternatively, one could simply run hydrody- 
namics event-by-event and compute the averages that are needed to compare to experiment 
[13, 15-19]. While event-by-event hydrodynamics is the best for this pragmatic purpose, the 
framework of non-linear response can yield valuable insight into the physics of these rather 
involved simulations. 



II. NON-LINEAR RESPONSE 
A. The cumulant expansion 

In hydrodynamic simulations of heavy ion collisions the medium is first modeled with 
an initial state Glauber model, then is evolved with hydrodynamics, and finally the particle 
spectrum is computed by making kinetic assumptions about the fluid. The final state particle 
spectrum for each event can be expanded in harmonics 



dN _ N 
dcb n 2ir 



n=l 



where here and below c.c. denotes complex conjugation. The root mean squares of v n are 
easily determined experimentally, and are given a special notation 

«»{2} = VW, (2-2) 

where ((. . .)) denotes the average over events. 

In the next sections we will describe how the momentum space response is related to the 
initial state geometry. To this end, the spatial distribution of the initial entropy density in 
the transverse plane, 

. , Tqs(x) 

is quantified with a cumulant expansion [9], where x = (x,y) = (r cos 0, r sin <ft) are the 
coordinates in the transverse plane and r G is the initial Bjorken time [20]. Specifically the 
n, m-th moment of the entropy distribution is defined as 

Pn , m = [ d 2 xp(x) ( r 2)(«-™)/2 r "v"W> , (2.4) 
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where (n — m)/2 is typically an integer. This moment is closely related to the n, m-th 
cumulant W„ )OT 

Wn,m Pn,m — contractions . (2.5) 

The meaning of eq. (2.5) will be clarified through examples, with additional details about 
the cumulant expansion relegated to the literature [9, 21]. The radial variation of p(x) is 
quantified by the radial cumulants, (r 2 ) and (r 4 ) — 2 (r 2 ) 2 , while the the azimuthal variation 
of p(x) is quantified by the azimuthal cumulants 



re 



e2e 



£3e 



«3<I>3 



(r 3 ) ' 
(r 3 ) 



(2.6) 
(2.7) 
(2.8) 



Here (. . .) denote an average over p(x) for a single event, and $i, $2 and $3 are the 
participant plane angles. These coordinate space angles are distinct from the momentum 
space angles ^i,^ 2 , and ^3. 

For the lowest harmonics the azimuthal cumulants and the azimuthal moments coincide, 
and these definitions will appear obvious to most readers. For the fourth harmonic and 
higher, we will depart from traditional moment based definition, and quantify the deforma- 
tions with cumulants rather than moments 2 



C 4 e 



i4<I>4 



(r^) - 3 (r 2 e^Y 



(2.9) 



The motivation for this definition can be seen by studying an elliptic Gaussian distribution, 

1 



p{x) 



2na x a y 



(2.10) 



which has C4 = 0, although €4 is non-zero and is of order e 2 . Similarly we define 

1 



(r 5 ) 



[(rV 5 *) - 10 (rV 2 *) (rV 3 ^)] 



(2.11) 



and remark that a Gaussian distribution deformed by an 63, 



S(X,T) OC 



1 + 



(r 3 ) e 3 
24 



d_ 

dx 



-3 



d Y d 



dy J dx 



g 2lT x 2tT H 



(2.12) 



has C5 = 0, although <(r 5 e* 5< ^ s ) is non-zero and of order 6263. 

We will characterize the hydrodynamic response to the cumulants defined above in the 
next section. 



2 For n>4we notate the cumulant based eccentricity by C n to differentiate this quantity from the moment 
based eccentricity e n . C n is equal to W n . n up to normalization and an overall factor of (r n ). 
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B. Non-linear response to the cumulants 



We expect the response of the system to be dominated by the lowest cumulants. Mo- 
tivated by Fourier analysis [9], we replace the general distribution p(x) with a Gaussian, 
eq. (2.10), whose second moments have been adjusted to reproduce (r 2 ) = a 2 x + cr 2 . and 
e 2 = (o"y — a x)/( a x + a y)- I n ^ ei '- we showed that a Gaussian + fourth order cumulants 
reproduces the results of smooth Glauber initial conditions in detail. If a Gaussian with a 
non-negligible €2 is simulated, the particle spectrum produced by this background contains 
all even harmonics 

™L = N_( 1 + W2 jW„-*2) + u, 4(ffl)e «<*.-*>> + . . . + c.c.) . (2.13) 

For small €2 the response coefficient W2 describes the linear response to the deformation and 
is proportional to e 2 , while ^4(22) describes the non-linear response and is proportional to 
e|. Below, we will assume that €2 is small enough that this scaling with €2 applies. Further, 
we have truncated the expansion in eq. (2.13) at quadratic order in €2, and will continue 
to do this implicitly from now on. The working assumption in this paper is that the most 
important non-linearity stems from the almond shape of the background. 

If the Gaussian distribution is perturbed by a small fourth order cumulant C 4 e* 4 * 4 , then 
the resulting particle spectra will be described by 

^ = ^ (1 + + W4 e^-^ + w m) e^-^ + c.c.) , (2.14) 

d(f)p 2n ' 

where captures the linear response to the fourth order cumulant and is proportional to 
C4 for small C4. In writing eq. (2.14) we have neglected terms proportional to £462, which 
can contribute to t>2 and reduce the perfect correlation between \l/ 2 and $ 2 - Comparing 
eq. (2.14) with the definition of V4, eq. (2.1), we see that v± is determined by the linear and 
quadratic response 

V4 e- 14 ^ = w A e- i4 ^ + w m) e- i4 ^ . (2.15) 
Squaring this result and averaging over events we see that 

M2} = {{vlf 2 = «Ke-^ + w m) e-^\ 2 f /2 . (2.16) 

In writing eq. (2.14) we have neglected the non-linear contributions of e\ and €3 to V4 since 
i>3 and vi are small compared to t; 2 for mid-peripheral collisions. 

Similarly, if the Gaussian background distribution is perturbed by a third order cumulant 
and a fifth order cumulant C5, then v$ is determined by a combination of the linear and 
non-linear response. The response to C5 is small [12], and therefore we will neglect the 
non-linearities due to €2^5, but we will keep the non-linearities due to €263. With this 
approximation scheme the particle spectrum through quadratic order reads 



dN _ N 
d(f) p 2tt 



w 1{ 2z)e i(t,p -^ +2 ^ + w 5(2 3)e i5 ^~ 3 * 3 " 2 * 2 + even harmonics + c.c) . (2.17) 
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Comparing this equation to the definition of V5, we see that 

v 5 {2} = ((\w 5 e- i5 ** + W5 ( 2 3)e-* (3 * 3+2 * 2) | 2 » 1/2 , (2.18) 

which is clearly analogous with t> 4 case. Finally, if the distribution has a net dipole asym- 
metry ei, then v\ is given a combination of the linear and non-linear response 

1 ,., , ^-i(3*3-2* 2 )|2\\ 1 /2 



Vl {2} = <(Ke— * + w im) e-^-^\Y , (2.19) 

where w± notates the linear response to e\. In writing this result for v\ we have neglected the 
non-linear interaction between V\ and V2, i.e. u>i(2i)- Thus eq. (2.19) makes the simplifying 
assumption that v\ is small compared to V3, while a more complete treatment would include 
a tyi(2i) contribution. 

Let us discuss how this formalism can be used to study the pr dependence of the flow. 
The particle spectra is expanded in harmonics 



dN dN 



dp T d(/)p dp T 



1 + Yl MPT)e in( ^-* niPT)) + c.c. , (2.20) 



n=l 



where the phase, *& n (pr), is i n general a function of pr- Then v n (pT){2} in the ^ n plane is 
normally defined as 

( ({»-(?r)"n cos(n(^n(pr)-^n)))) « > 1 

Vn{PT){2} = { ((vi(pt)vi co"s(^i(p T )-^i))) _ n > (2-21) 

vi{2} 



n = 1 



where we have inserted an extra minus sign for t>i(Pr)> since the integrated V\ is negative. 
The phase angle ^ n {PT) is often assumed to equal \l/ n . Using the formalism outlined above 
we write Vi{px) as a sum of the linear and non-linear response 

MPT)e-** l(PT) = w.Me-^ + w 1(23) (p T )e- i3 ^ +i2 ^ . (2.22) 

Then the numerator of V\(pt){2} is given by 

((vi(Pt)vi cos(#i(p T ) - ^i)» = 

I^iWibi + wi(23)(Pr)wi(23) + [«'i(Pr)«'i(23) + ^1(23) (pr)wi] cos($i - 3$ 3 + 2$ 2 )^) , 

(2.23) 

and the denominator is given by the integrated expression for t>i{2}, eq. (2.19). Similar 
expressions follow for Uj(pt) and 1*5 (pr)- 

Finally, let us place some older measurements and calculations of v<±(pt) into context 
[22-28]. Traditionally, what was referred to as v^pr) would today be called v^ipr) in the 
^2 plane: 

/ \ ro1 _ (Mpr)^cos(4fr 4 (p T ) - 2^ 2 - 2tt 2 ))) 
V4(22)[Pt){2\ = — j-r • (2.24) 

^21^J 

As discussed in the conclusions, the differences between v ^22) (pt) {2} and v 4 (p T ){2} can be 
used to partially disentangle the linear and non-linear response. 
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C. Summary 



The goal of the present work is to compute the linear and non-linear response coefficients, 
and to use these coefficients together with an initial Glauber model to determine {{v^}) with 
Eqs. (2. 16), (2. 18), and (2.19). For v$ the step by step procedure is: (i) use hydrodynamics 
to determine the response coefficients 

f, a^d (2.25) 

for vanishingly small C5 and €263; (ii) use a Glauber model to determine the geometric coef- 
ficients that are needed in eq. (2.18), ((Cf)), (((e 2 e 3 ) 2 )), and ((C 5 e 2 e 3 cos(5$ 5 - 3$ 3 - 2$ 2 ))); 
(iii) combine these results in eq. (2.18) to determine the complete hydrodynamic prediction 
for ((vl)). The necessary Glauber correlations are determined using the Phobos Monte Carlo 
Glauber Model [29], and we note that there is a very strong geometric correlation between 
participant planes differing by two, e.g. 

((Cr e 2 e 3 e^- 3 ^- 2 ^)) , and ((exW'"--*-^)) . (2.26) 

This geometric correlation can be studied analytically in an independent source model [30] , 
and is easily attributed to the elliptic shape of the overlap region [9, 30, 31]. 



III. HYDRODYNAMIC SIMULATIONS 
A. Ideal and viscous hydrodynamics 

To calculate the non-linear response we use a hydrodynamics code that implements con- 
formal second order hydrodynamics [32] . The numerical scheme is based on a central scheme 
developed and tested in Ref. [33], although the equations of motion for the are some- 
what different from what was studied in that work 3 , rj/s is held constant, and the ra- 
tio of second order hydro parameters are taken from their AdS/CFT values [32, 34], e.g. 
T^/ijri/sT) = 4 — 2 In 2. The equation of state partially parametrizes lattice results and 
was used previously by Romatschke and Luzum [35]. Finally, we have followed the time 
"honored" constant temperature freezeout prescription, with Tf D = 150 MeV. For simplicity 
we have adopted the popular quadratic ansatz for the viscous correction to the thermal 
distribution function [2] 

f(P) = UP) + 6f(P) , 5f(P) = 2fe\%)°T^ P ^ ' (3J) 

where f (P) = l/(exp(— P ■ U(X)/T) =F 1) is the equilibrium distribution, e + V is the 
enthalpy, and 5f is the first viscous correction [2, 36]. Although we have used the quadratic 
ansatz in this work, a linear ansatz is probably more appropriate for QCD-like theories and 
can effect the integrated flow for the higher harmonics [25, 37]. 

3 However, when additional non-conformal second order gradients are added to our equations of motion 
and the parameters are matched, our current numerical can be compared directly to Ref. [33]. If this is 
done, the two hydro-codes yield the same answers to 0.1% for the type of problems considered in this 
work. 
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For the simulations shown below we have followed the centrality classification given in 
Ref. [13] which is documented in Table I. of that work. Our procedure to determine the 
response coefficient at a given impact parameter largely follows Ref. [9], which should be 
referred to for additional details - see especially Appendix A of that work. Briefly, for each 
impact parameter we determine the average squared radius (r 2 ), and initialize a Gaussian 
distribution that is deformed by the appropriate cumulant. The Gaussian is normalized to 
reproduce the total entropy in the event. For instance, to determine the w^^3) we initialize 
the distribution given in eq. (2.12) with 62 = £3 = 0.02. A technical complication is that the 
distribution in eq. (2.12) must be regulated [9], and the regularization procedure introduces a 
small C5. However, the spurious C5 decreases faster than e 3 and can be made arbitrarily small 
compared to the signal. Empirically we find that the spurious C5 decreases approximately 
as e 5 , and the v§ from the spurious cumulant is negligibly small compared to the v$ from 
the £263 combination. 



B. The non-linear response coefficients in ideal and viscous hydrodynamics 

In this section we will study the non-linear response coefficients systematically. In par- 
ticular we study how the linear and non-linear response coefficients depend (i) transverse 
momentum, (ii) centrality, and (iii) shear viscosity. 



1. Momentum dependence of the response coefficients 

Fig. 1 examines the pt dependence of the linear and non-linear response coefficients, W4 
and u>4(22)> which are characteristic of the response coefficients more generally. First, focus 
on the ideal curves in Fig. 1(a) and (b). At large pt the non-linear response curves show a 
characteristic quadratic rise with px, while the linear response curves show a characteristic 
linear rise. This difference between the non-linear and linear response is known from previous 
studies of [23]. Later, when examining non-linear corrections to v\ (see Fig. 6), we will see 
that the non-linear corrections are most important at high px and exhibit a characteristic 
quadratic rise. Comparing Fig. 1(a) and (b), we see that viscous corrections are smaller 
for the non-linear response ^4(22)(PT)/ e 2' than for the linear response iy 4 (p T )/C 4 . This is a 
generic result as will be discussed in detail in Section IIIB3. 

We also note that the linear response curves shown in Fig. 1(a) change sign for sufficiently 
large viscosity. This is an artifact of the first viscous correction, 5f, and the quadratic ansatz. 
To see this, we have plotted w±{px) and and w^22){pt) using only the unmodified distribution 
function f in Fig. 1(a) and (b). For large viscosity the 5f correction to V4 and v 5 is not small 
compared to the ideal contribution f , and this causes a reduction of the response, which 
is more pronounced for the higher harmonics, t>4 and v$. In full kinetic theory calculations 
W4/C4 and W5/C5 remain positive and approach zero as the viscosity is increased [12]. Thus, 
the negative u> 4 /C 4 indicates that the first viscous correction has become too large to be 
trusted. Below, we will simply set the response coefficients to zero when this is the case. 
Experience with kinetic theory suggests that this ad hoc procedure is not far from what 
really happens. 
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FIG. 1. (color online) The linear and non-linear response coefficients for w^{pt) and 104(22) (j>t), 
for ideal and viscous hydrodynamics. The grey curves (which are shown only for r\j s = l/4ir) ex- 
hibit the resulting response when the equilibrium distribution f Q is used, and the viscous correction, 
5f, is neglected (see eq. (3.1)). 



2. Centrality dependence of the response coefficients 



Fig. 2 shows the linear and non-linear response coefficients in ideal and viscous hydro- 
dynamics. There are several salient features contained in these plots. First, note that the 
magnitude of the linear response coefficient W5/C5 is quite small in the viscous case, and 
W5/C5 has been multiplied by ten to make the curves visible. The non-linear response w^23) 
coefficient is significantly larger. The implications of this difference will be studied in the 
next section when we multiply the response coefficients by C 5 and e 2 e 3 respectively. Second, 
all of the response coefficients are reduced by viscosity, especially in non-central collisions. 

The viscous and curves stop abruptly as a function of centrality, since we 

have truncated the curves when response falls below zero. As discussed above (see Fig. 1), 
this is because viscous corrections to the thermal distribution function (Sf) become larger for 
more peripheral collisions, and this correction is magnified by the high harmonic number. We 
have therefore truncated the and w§ response curves when the response turns negative. 
At this point 5f constitutes an order one correction and can no longer be trusted. 



3. Dependence on viscosity 



It is interesting to note that viscous reduction for Wi/ei is smaller than for w 4 /C 4 and 
W5/C5. This pattern of viscous corrections for linearized perturbations is studied further 
in Fig. 3(a). Each linearized perturbation labeled by n, m-th cumulant is damped by a 
factor ~ exp(— r n m Tfi na i) relative to ideal hydrodynamics, where Tfi na i is an estimate for the 
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FIG. 2. The linear and non-linear response coefficients for ideal and viscous hydro. In the viscous 
case the curves are truncated when the response coefficients turn negative, i.e. outside of the 
regime of validity of viscous hydro. 



duration of the event. Analytical work shows that the damping coefficients r„ jjra scale as 



f)i,m Tfinal 



mfp / H fTl 



L 



+ m 



(3.2) 



for a conformal equation of state and a particular background flow [21]. Thus, each power of 
r 2 and each harmonic order in eq. (2.4) increases (n—m)/2+m by one unit. Our numerical 
work (Fig. 3(a)) is not limited to the conformal equation of state or the particular background 
flow of Ref. [21], and shows that this scaling is reasonably generic [12, 38]. Specifically, the 
formal estimate given in eq. (3.2) implies a definite pattern among the viscous corrections 
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FIG. 3. (a) Linear response coefficients function of viscosity relative to the ideal response, 

(b) A comparison of non-linear and linear response coefficients as a function of viscosity, e.g. 1^5(23) 
records the V5 produces by a combination of e 2 and £3. The negative values for large viscosity are 
spurious, and lie beyond the region of applicability of viscous hydrodynamics. 



to v n : 

Awi Aw 2 .77 Aw 3 77 Aw A rj Aw 5 77 

rr- ~ 7-r- OC 4 - , rr- OC 9 - , OC 16 - , rr- OC 25 - . 

Wf Wf S Wf S Wf S wf S 

(3.3) 

where Aw = w vlscous — u> ldeal , and w ld is the ideal hydro response coefficient. Note, in 
particular that the viscous corrections v i and v 2 are similar since V\ and v 2 respond to the 
dipole asymmetry, W^i, and the ellipticity, 1^2,2, respectively [38]. Since the slopes of the 
V\ '■ v 2 : t>3 : V4 : v 5 curves in Fig. 3(a) have approximately the expected ratios 4 : 4 : 9 : 16 : 
25, our numerical work qualitatively confirms this pattern of viscous corrections. 

Fig. 3(b) compares the damping rate for the non-linear response coefficients to the corre- 
sponding linear response coefficients. Take ^5(23) for example. Since ^5(23) is of order t> 2 t>3 
we expect the damping of this non-linear perturbation to scale as ~ e~ r ' 2 - 2T e~ r3 ' 3T , and thus 
the damping rate r 5 ( 2 3) is expected to scale as 

T5(23) ~ r 2i 2 + r 3i3 . (3.4) 

Thus, we expect the non-linear and linear response coefficients for v$ to scale as 

-^MoclS^, -^oc25^. (3.5) 

W 5(23) S W f S 

Comparing the slopes of the non-linear and linear response curves in Fig. 3(b), we see that 
the slope of the Aw 5 ^)/w^ 2 ^ curve is approximately half of the corresponding Aw 5 /wf, 
and is qualitatively consistent with our heuristic estimate of 13/25. w 4 ( 2 2) and show a 
similar pattern of viscous corrections. Finally our estimates seem only partially applicable to 
V\. For instance, the reasoning of eq. (3.4) predicts that the non-linear damping rates, Tn2z) 
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and r 5 (23), should be equal. However, the slope of Awi( 2 3) / wl U23) * s significantly smaller than 
the Aw 5 (23)/w5 t ( 2 3- ) , and contradicts this reasoning. Clearly, the non-linear viscous damping 
of Vi is a special case which will have to be investigated more completely at a later date. 

IV. RESULTS AND DISCUSSION 
A. Results 

Having clarified the non-linear hydrodynamic response, we study the phenomenological 
implications of these response coefficients. Fig. 4 shows v\, V4 and vc, including the linear 
and non-linear response as outlined in Section II, and is the principal result of this work. 

Examining this figures we see that the non-linear response is an important correction 
for v\, and essential for V4 and V5. The contribution of the non-linear response to the 
total flow increases towards peripheral collisions, and for t> 4 and v§ is of order 50% in mid- 
peripheral collisions. This is roughly compatible with simulation results from event-by-event 
hydrodynamics [13, 14]. Especially for viscous hydrodynamics and for V5, the linear response 
is negligible in all but the most central bin. Even in the most central bin, the non-linear 
contribution to v 5 is about 50% of the total. It is notable, if expected, that for v\ viscosity 
reduces the non-linear contribution relative to the total, while for V4 and v§ viscosity increases 
the non-linear contributions. This is consistent with the discussion given in Section II B. 

It will be quite interesting to measure the complete set of event planes ^3, ^4, ^5) 

and their inter-correlations. These measurements will place a strong experimental constraint 
on the relative of importance of the non-linear response [11]. For example, if the non-linear 
response is dominant (as implied by the viscous v§ curves), then a stronger than geometric 
correlation is expected for certain experimental averages, e.g. (cos (5^5 — 3^/ 3 — 2^/ 2 ))- 

Next we examine the pr dependence of the t>i,t>4 and v$. Since U4 and v$ are dominated 
by the non-linear response we will present our results by scaling v\ and t>2t>3 respectively. 
Many of the points raised in this and the next paragraph are familiar from earlier studies 
of t>4 in the ^2 plane. In particular, the importance of non-linearities and fluctuations in 
determining the experimental v±jv\ ratio was understood previously [23, 24]. 

First we note that according to an old argument by Borghini and Ollitrault [23], v^jv\ 
should approach 1/2 at large momentum in ideal hydrodynamics for any given event due to 
the non-linearities inherent in the phase space distribution. Their result is easily generalized 
to t>5, t>5 = t>2^3- The argument follows by computing the freezeout distribution in a saddle 
point approximation [39], and can be schematically understood by examining the thermal 
factor in an approximately radially symmetric flow profile. The transverse flow vector as a 
function of the spatial azimuthal angle relative to the reaction plane is 

ut = (u x ,u y ) ~ (ut(4>) cos0, ut{4>) sin0) , (4.1) 

where in the second step we have assumed that the flow is approximately radially symmetric. 
The transverse flow velocity is then expanded in harmonics 

ut{4>) = + 2u£ cos 2(f) + cos 40 + other harmonics , (4.2) 
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FIG. 4. vi, v& and v§ versus centrality in ideal and viscous hydrodynamics. To keep the ideal and 
viscous curves on the same scale we have multiplied the viscous V5 curves by a factor of two. In 
the viscous case, the linear response is neglected when the response coefficients turn negative, i.e. 
outside of the region of applicability of viscous hydrodynamics. 
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and the thermal factor with p = {j>t cos p , pr sin P ) reads 



-'(cos 4> p —<f>) 



1 , 2 PT (2) 

i + -—uy cos. 

2p 



( 2 ) nj, 
—My COS 20 



+ 



1 , ^ T (2) , /PT (2)V AM 

1 + ~jT U T cos20 p + ^— u T ; J cos40 p + 



(4.3) 
(4.4) 



The leading exponential strongly correlates coordinate space angle and the momentum 
space angle <f> p . In the second line we have anticipated the saddle point approximation, 
(which realizes this correlation) and set (f) ~ <j) p in the post-exponent. The second term in 
square brackets determines the linear response coefficient w 2 and rises linearly with momen- 

(2) 

turn, w 2 ~ ptu t /T. The third term determines the non-linear response coefficient u>4(22); 
and grows quadratically with momentum, u>4(22) ~ \{Ptu^ /T) 2 . At high pt this quadratic 
growth overwhelms the (neglected) linear response due to v!^\ and leads to the characteristic 



relation v± 



2 u 2 . 



An entirely identical argument shows that v 5 = v 2 v 3 at high momentum 



in ideal hydrodynamics. 

The Borghini-Ollitrault argument given above shows that the response coefficients in 
ideal hydro should asymptote at large momentum, 



^4(22)/^ 

(w 2 /e 2 ) 2 



w^ 23) /{e 2 e 3 ) 
{w 2 /e 2 ){w 3 /e 3 ) 



-> 1 . 



(4.5) 



When fluctuations are included these asymptotic relations are modified [24]: 



M2} 
v 2 {2} 2 



Pt- 



1 

>oo 2 



4) 



(4Y 



1/2 



v 5 {2} 



v 2 {2}v 3 {2} 



Pt- 



((^26 3 ) 2 ) 

(4) (el) 



1/2 



(4.6) 



Previous studies of V4 in the ty 2 plane (see Section II B) have shown that such geometrical 
factors are essential to reproducing the centrality dependence of v±jv\ [24]. The following 
table records the geometrical ratios in eq. (4.6) as a function of centrality. 



Centrality % 


2.5 


7.5 


12.5 


17.5 


25.0 


35.0 


45.0 


55.5 


\l(4)/(4? 


1.40 


1.33 


1.26 


1.22 


1.20 


1.18 


1.17 


1.16 


V / ((^3) 2 )/((ei) (el)) 


0.99 


0.97 


0.96 


0.95 


0.94 


0.93 


0.930 


0.92 


V(4)/(4) 


2.12 


2.78 


3.22 


3.46 


3.56 


3.40 


3.04 


2.64 



We have found that rather large px is needed to see the non-linear limit given by eq. (4.6). 
In the current framework, the linear and non-linear response terms, and their interference, 
determine the full result 



((w 4 (p T )w 4 + U>4(22)(PrV 4 (22) + [w 4 (p T ) ^4(22) + W 4(22) (p T ) W 4 ] COs4($ 4 - $ 2 ))) 

M2|(pt) " M2} 

(4-7) 

Fig. 5 shows the complete result for v 4 {2} / v 2 {2} 2 (scaled by J (e|) / (e 2 .) 2 ) for ideal and 
viscous hydrodynamics. Focusing on the ideal results, we see that full results (the solid 
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FIG. 5. Results for V4 and v§ for ideal and viscous hydrodynamics at various impact parameters. 
The Borghini-Ollitrault expectation is indicated by the arrows for the ideal V4 and v§ curves [23]. 



lines) approach the non-linear expectation of Borghini and Ollitrault (the dashed line) only 
very slowly. This is in large part because w^(pt) is only qualitatively linear at sub-asymptotic 
Pt and increases almost quadratically at intermediate pt ~ 1.5 GeV, momentarily keeping 
up with the non-linear response. When viscous corrections are included, the non-linear 
results become dominant in peripheral collisions. Similar results for t> 5 in ideal and viscous 
hydrodynamics are also shown in Fig. 5. In the viscous case, the non-linear result gives 
almost the full t>5{2} for all centrality classes shown. 

It is worth noting that the magnitude of the viscous corrections as a function of pt for 
t>4 and V5 are sensitive to ansatz used for the viscous distribution function, 5f [26]. In 
particular, the quadratic ansatz used in this work assumes that the quasi-particle energy 
loss is independent of momentum, dp/dt oc const. A linear ansatz for 5f is better motivated 
for QCD like theories and results in smaller viscous corrections for w 4 and v§ as a function 
of pt [37]. A complete discussion of this point is reserved for future work. 
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FIG. 6. vi(pt){2} (eq. (2.21)) in ideal and viscous hydrodynamics from the linear response to e±, 
the non-linear response to €263, and the total response, eq. (2.23). 

Fig. 6 presents the corresponding analysis for Vi(pt). We see that the non-linear terms 
provide a correction to the linear response which grows with pt due to the quadratic de- 
pendence of the non-linear response coefficients, u>i(23) oc p\. We note that the viscous 
corrections are approximately the same for t>i(pr) and V2(pr), as expected from the discus- 
sion of viscous corrections given in Section II B. 

B. Discussion 

We have presented a framework of non-linear response to understand the higher harmonics 
generated in heavy ion collisions. Then we extracted the non-linear response coefficients 
using ideal and viscous hydrodynamics and studied the dependence on the shear viscosity, 
in Fig. 2. The pattern of viscous corrections is further analyzed in Fig. 3 and explained 
in Section III B. Generally, when the harmonic order is large, the non-linear response is 
less damped than the corresponding linear response. Thus, when viscosity is included in 
hydrodynamic simulations, the non-linear response becomes increasingly important for the 
higher harmonics. This qualitative reasoning is confirmed in Fig. 4 which shows v\, V4 and 
v§ using linear and non-linear response and is the principal result of this work. We see that 
the non-linear response is essential for V4 and v§, and constitutes an important correction 
for V\. 

Experimentally, the relative contributions of the linear and non-linear response can be 
disentangled by measuring v§ in the 2^/2 + 3\&3 and planes, i.e. by measuring 

^5(23) = (cos(50 p - 2^ 2 - 3*3)) , and v 5(5) = (cos5(0 p - * 5 )> • (4.8) 

Although a full discussion of this and similar measurements is reserved for future work, 
a qualitative expectation based on Fig. 4(e) and (f) is that the (cos(5\I / 5 — 3^% — 2^2)) 
correlation should be strong compared to the geometric average, and should change rapidly 
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from central to mid-central collisions. Qualitatively, this is precisely what was observed 
recently by the ATLAS collaboration [11]. 

The non-linear response can also be studied by analyzing the pt dependence of the flow 
harmonics. Fig. 5 and Fig. 6 exhibit v±(pt), v$(pt), and vi(pr)- In ideal hydrodynamics at 
large px we expect to find V4 = \v\ on an event by event basis. Our non-linear response 
coefficients corroborate this non-linear expectation for t> 4 and an analogous relation for t> 5 , 
V5 = V2V3. However, since what is normally measured is f4{2}/(u 2 {2}) 2 and not ((v^/vf)), 
this ideal non-linear expectation must be multiplied by ((e^) / (c^) 2 ) 1 ^ when comparing to 
the experimental data [24]. In addition, this expectation of ideal hydrodynamics is broken 
by viscous corrections, and by the linear response to the fourth order cumulant C4 {i.e. €4). 
When all of these corrections are taken into account, we find that relations such as V4 = \v\ 
and t> 5 = v 2 v 3 provide only a rough guide to the full result. 

Throughout we have assumed perfect correlation between \l/ 2 and $2 and ^3 and $3. 
This strict correlation is only approximately true. For instance the combination of a v\ and 
a v 3 can yield a v 2 , 

V2 e~ i2 ^ = w 2 e~^ + w 2{13) e-^ +l ^ . (4.9) 

This naturally provides a correlation between the ^2 arid ^3 plane, although the geometric 
correlation between $2 and $3 is negligibly small. Indeed the (^2, ^3) correlation, which was 
very recently observed by the ATLAS collaboration [11], is too large to be easily explained 
with the geometric correlations of the Glauber model. Similarly, assuming that the linear 
response to e 6 is negligible, one could expect that in central collisions w 6 is determined by the 
quadratic response to v 3 , while in peripheral collisions vq is determined by a cubic response 
to v 2 

v s e- im * = w^e^ 2 + u^e"** 3 . (4.10) 

Qualitatively, this pattern is consistent with the observed (\l/ 6 , ^/ 3 ) and (\l/ 6 , \& 2 ) correlations 
presented in [11]. It will be interesting to see if all of the observed correlations can be 
quantitatively understood with the non-linear response theory outlined in this paper. A full 
quantitative comparison with the experimental data is reserved for future work. 
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